
## @knitr MODELS_SETUP_SUBSTANTIVE_EFFECTS

rm(list=ls())

#library(foreign)
#library(stringr)
#library(xlsx)

# from \url{https://www.r-bloggers.com/identifying-the-os-from-r/}
get_os <- function(){
  sysinf <- Sys.info()
  if (!is.null(sysinf)){
    os <- sysinf['sysname']
    if (os == 'Darwin')
      os <- "osx"
  } else { ## mystery machine
    os <- .Platform$OS.type
    if (grepl("^darwin", R.version$os))
      os <- "osx"
    if (grepl("linux-gnu", R.version$os))
      os <- "linux"
  }
  tolower(os)
}

get_machine  <- function(){
    sysinf  <- Sys.info()
    if (!is.null(sysinf)){
        my.machine  <- sysinf['nodename']
    }
    tolower(my.machine)
}

if(get_os()=="windows" & get_machine()=="mompracen") {
    path <- setwd('C:/Users/giaco/Documents/MEGAsync/')
} else if (get_os()=="windows" & get_machine()!="mompracen") {
    path <- setwd('C:/Users/grieco/Dropbox/')
} else if (get_machine()=="aramis") {
    path <- setwd('/home/giacomo68/Dropbox/')
} else {
    path <- setwd('/home/giacomo/Dropbox/')
}


library(survival)
## library(coxme) 
library(xtable)
library(ggplot2)
library(english)
library(Publish)
library(showtext)
library(extrafont)
library(mice)
library(mitools)
library(sandwich)
library(lmtest)

Data <- read.csv('gc-jg-project/kampala/data/preliminary-data-kampala.csv',
                 stringsAsFactors=FALSE)


Data$country_id <- as.numeric(factor(Data$country_name))

m <- 10

imp <- mice(Data[,c("kampala",
                 "sidea.1991",
                 "sideb.1991",
                 "my.mil.cap.lag",
                 "v2x_polyarchy.lag",
                 "v2x_libdem.lag",     
                 "v2x_partipdem.lag",
                 "v2x_delibdem.lag", 
                 "v2x_egaldem.lag",
                 "gdp.cap.lag",
                 "pop.lag",
                 "common.law",
                 "year",
                 "rule.law.lag",
                 "aid.recipient",
                 "aid.recipient.q",
                 "All_NGO.lag",
                 "country_id")], m = m, maxit = 10, seed = 1510)

models <- list()
model_results <- list()
merged_imputations <- list()
m1.model_results <- list()
m2.model_results <- list()
m3.model_results <- list()
m4.model_results <- list()
m5.model_results <- list()
m6.model_results <- list()

for (i in 1:m) {
    complete_data <- complete(imp, i)

    complete_data$.imp <- i
    
    merged_data <- merge(complete_data, 
                         Data[,c("country_id","year","country_name",
                                 "kampala.t0","kampala.t1",
                                 "sidea.1945",
                                 "sideb.1945",
                                 "sidea.force.1945",
                                 "sideb.force.1945",
                                 "sidea.no.force.1945",
                                 "sideb.no.force.1945",
                                 "sidea.force.1991",
                                 "sideb.force.1991",
                                 "sidea.no.force.1991",
                                 "sideb.no.force.1991"
                                 )], 
                        by = c("country_id", "year"),  # Note: quoted variable names
                        all.x = TRUE)  # Keep all rows from complete_data
    
    merged_imputations[[i]] <- merged_data


    current_data <- merged_imputations[[i]]

    m1 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
              sidea.1991
            + sideb.1991
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_id)
           ,data=current_data)

     m1.vcov_cluster <- vcovCL(m1, cluster = current_data$country_id)

     m1.model_results[[i]] <- list(
        coefficients = coef(m1),
        variance = m1.vcov_cluster,
        n = nrow(current_data)
       )
    

    
    m2 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                sidea.1945
            + sideb.1945
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_name)
           ,data=current_data)

    m2.vcov_cluster <- vcovCL(m2, cluster = current_data$country_id)

    m2.model_results[[i]] <- list(
        coefficients = coef(m2),
        variance = m2.vcov_cluster,
        n = nrow(current_data)
       )
    


    m3 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.force.1991
                + sideb.force.1991
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)

    m3.vcov_cluster <- vcovCL(m3, cluster = current_data$country_id)

    m3.model_results[[i]] <- list(
        coefficients = coef(m3),
        variance = m3.vcov_cluster,
        n = nrow(current_data)
       )
    
    

    m4 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.force.1945
            + sideb.force.1945
            + my.mil.cap.lag
            + v2x_libdem.lag
            + log(gdp.cap.lag)
            + log(pop.lag)
            + common.law
            + cluster(country_name)
           ,data=current_data)

    m4.vcov_cluster <- vcovCL(m4, cluster = current_data$country_id)

    m4.model_results[[i]] <- list(
        coefficients = coef(m4),
        variance = m4.vcov_cluster,
        n = nrow(current_data)
       )
    


    m5 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.no.force.1991
                + sideb.no.force.1991
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)
    
    m5.vcov_cluster <- vcovCL(m5, cluster = current_data$country_id)

     m5.model_results[[i]] <- list(
        coefficients = coef(m5),
        variance = m5.vcov_cluster,
        n = nrow(current_data)
       )
    


    m6 <- coxph(Surv(kampala.t0,kampala.t1,kampala)~
                    sidea.no.force.1945
                + sideb.no.force.1945
                + my.mil.cap.lag
                + v2x_libdem.lag
                + log(gdp.cap.lag)
                + log(pop.lag)
                + common.law
                + cluster(country_name)
               ,data=current_data)
    
    m6.vcov_cluster <- vcovCL(m6, cluster = current_data$country_id)

    m6.model_results[[i]] <- list(
        coefficients = coef(m6),
        variance = m6.vcov_cluster,
        n = nrow(current_data)
       )
    
}


m1.pooled_results <- MIcombine(
    results = lapply(m1.model_results, function(x) x$coefficients),
    variances = lapply(m1.model_results, function(x) x$variance)
)

m2.pooled_results <- MIcombine(
    results = lapply(m2.model_results, function(x) x$coefficients),
    variances = lapply(m2.model_results, function(x) x$variance)
)

m3.pooled_results <- MIcombine(
    results = lapply(m3.model_results, function(x) x$coefficients),
    variances = lapply(m3.model_results, function(x) x$variance)
)

m4.pooled_results <- MIcombine(
    results = lapply(m4.model_results, function(x) x$coefficients),
    variances = lapply(m4.model_results, function(x) x$variance)
)

m5.pooled_results <- MIcombine(
    results = lapply(m5.model_results, function(x) x$coefficients),
    variances = lapply(m5.model_results, function(x) x$variance)
)

m6.pooled_results <- MIcombine(
    results = lapply(m6.model_results, function(x) x$coefficients),
    variances = lapply(m6.model_results, function(x) x$variance)
)



# View final results
summary(m1.pooled_results)
summary(m2.pooled_results)
summary(m3.pooled_results)
summary(m4.pooled_results)
summary(m5.pooled_results)
summary(m6.pooled_results)


coef.m1 <- m1.pooled_results$coefficients
se.m1 <- sqrt(diag(m1.pooled_results$variance))
z.m1 <- coef.m1/se.m1
p.m1 <- 2*pnorm(-abs(z.m1))
conf.m1 <- confint(m1.pooled_results)
m.m1 <- cbind(coef.m1,se.m1,p.m1,conf.m1)

coef.m2 <- m2.pooled_results$coefficients
se.m2 <- sqrt(diag(m2.pooled_results$variance))
z.m2 <- coef.m2/se.m2
p.m2 <- 2*pnorm(-abs(z.m2))
conf.m2 <- confint(m2.pooled_results)
m.m2 <- cbind(coef.m2,se.m2,p.m2,conf.m2)

coef.m3 <- m3.pooled_results$coefficients
se.m3 <- sqrt(diag(m3.pooled_results$variance))
z.m3 <- coef.m3/se.m3
p.m3 <- 2*pnorm(-abs(z.m3))
conf.m3 <- confint(m3.pooled_results)
m.m3 <- cbind(coef.m3,se.m3,p.m3,conf.m3)

coef.m4 <- m4.pooled_results$coefficients
se.m4 <- sqrt(diag(m4.pooled_results$variance))
z.m4 <- coef.m4/se.m4
p.m4 <- 2*pnorm(-abs(z.m4))
conf.m4 <- confint(m4.pooled_results)
m.m4 <- cbind(coef.m4,se.m4,p.m4,conf.m4)

coef.m5 <- m5.pooled_results$coefficients
se.m5 <- sqrt(diag(m5.pooled_results$variance))
z.m5 <- coef.m5/se.m5
p.m5 <- 2*pnorm(-abs(z.m5))
conf.m5 <- confint(m5.pooled_results)
m.m5 <- cbind(coef.m5,se.m5,p.m5,conf.m5)

coef.m6 <- m6.pooled_results$coefficients
se.m6 <- sqrt(diag(m6.pooled_results$variance))
z.m6 <- coef.m6/se.m6
p.m6 <- 2*pnorm(-abs(z.m6))
conf.m6 <- confint(m6.pooled_results)
m.m6 <- cbind(coef.m6,se.m6,p.m6,conf.m6)


#### Substantive effects

pi.m1.mida.10 <- round(100/(1+exp(coef.m1[1]*(1-0))),2)
ci.m1.mida.10 <- round(100/(1+exp(conf.m1[1,]*(1-0))),2)
sub.m1.mida.10 <- c(pi.m1.mida.10,ci.m1.mida.10[2],ci.m1.mida.10[1])

pi.m1.mida.20 <- round(100/(1+exp(coef.m1[1]*(2-0))),2)
ci.m1.mida.20 <- round(100/(1+exp(conf.m1[1,]*(2-0))),2)
sub.m1.mida.20 <- c(pi.m1.mida.20,ci.m1.mida.20[2],ci.m1.mida.20[1])

pi.m1.mida.30 <- round(100/(1+exp(coef.m1[1]*(3-0))),2)
ci.m1.mida.30 <- round(100/(1+exp(conf.m1[1,]*(3-0))),2)
sub.m1.mida.30 <- c(pi.m1.mida.30,ci.m1.mida.30[2],ci.m1.mida.30[1])

pi.m1.midb.10 <- round(100/(1+exp(coef.m1[2]*(1-0))),2)
ci.m1.midb.10 <- round(100/(1+exp(conf.m1[2,]*(1-0))),2)
sub.m1.midb.10 <- c(pi.m1.midb.10,ci.m1.midb.10[2],ci.m1.midb.10[1])

pi.m1.midb.20 <- round(100/(1+exp(coef.m1[2]*(2-0))),2)
ci.m1.midb.20 <- round(100/(1+exp(conf.m1[2,]*(2-0))),2)
sub.m1.midb.20 <- c(pi.m1.midb.20,ci.m1.midb.20[2],ci.m1.midb.20[1])

pi.m1.midb.30 <- round(100/(1+exp(coef.m1[2]*(3-0))),2)
ci.m1.midb.30 <- round(100/(1+exp(conf.m1[2,]*(3-0))),2)
sub.m1.midb.30 <- c(pi.m1.midb.30,ci.m1.midb.30[2],ci.m1.midb.30[1])

pi.m1.poly.24 <- round(100/(1+exp(coef.m1[4]*(4-2))),2)
ci.m1.poly.24 <- round(100/(1+exp(conf.m1[4,]*(4-2))),2)
sub.m1.poly.24 <- c(pi.m1.poly.24,ci.m1.poly.24[2],ci.m1.poly.24[1])

pi.m1.poly.26 <- round(100/(1+exp(coef.m1[4]*(6-2))),2)
ci.m1.poly.26 <- round(100/(1+exp(conf.m1[4,]*(6-2))),2)
sub.m1.poly.26 <- c(pi.m1.poly.26,ci.m1.poly.26[2],ci.m1.poly.26[1])

pi.m1.poly.28 <- round(100/(1+exp(coef.m1[4]*(8-2))),2)
ci.m1.poly.28 <- round(100/(1+exp(conf.m1[4,]*(8-2))),2)
sub.m1.poly.28 <- c(pi.m1.poly.28,ci.m1.poly.28[2],ci.m1.poly.28[1])

pi.m1.gdp.24 <- round(100/(1+exp(coef.m1[5]*(log(4000)-log(2000)))),2)
ci.m1.gdp.24 <- round(100/(1+exp(conf.m1[5,]*(log(4000)-log(2000)))),2)
sub.m1.gdp.24 <- c(pi.m1.gdp.24,ci.m1.gdp.24[2],ci.m1.gdp.24[1])

pi.m1.gdp.26 <- round(100/(1+exp(coef.m1[5]*(log(8000)-log(2000)))),2)
ci.m1.gdp.26 <- round(100/(1+exp(conf.m1[5,]*(log(8000)-log(2000)))),2)
sub.m1.gdp.26 <- c(pi.m1.gdp.26,ci.m1.gdp.26[2],ci.m1.gdp.26[1])

pi.m1.gdp.212 <- round(100/(1+exp(coef.m1[5]*(log(12000)-log(2000)))),2)
ci.m1.gdp.212 <- round(100/(1+exp(conf.m1[5,]*(log(12000)-log(2000)))),2)
sub.m1.gdp.212 <- c(pi.m1.gdp.212,ci.m1.gdp.212[2],ci.m1.gdp.212[1])

pi.m1.pop.410 <- round(100/(1+exp(coef.m1[6]*(log(10000)-log(4000)))),2)
ci.m1.pop.410 <- round(100/(1+exp(conf.m1[6,]*(log(10000)-log(4000)))),2)
sub.m1.pop.410 <- c(pi.m1.pop.410,ci.m1.pop.410[2],ci.m1.pop.410[1])

pi.m1.pop.420 <- round(100/(1+exp(coef.m1[6]*(log(20000)-log(4000)))),2)
ci.m1.pop.420 <- round(100/(1+exp(conf.m1[6,]*(log(20000)-log(4000)))),2)
sub.m1.pop.420 <- c(pi.m1.pop.420,ci.m1.pop.420[2],ci.m1.pop.420[1])

pi.m1.pop.450 <- round(100/(1+exp(coef.m1[6]*(log(50000)-log(4000)))),2)
ci.m1.pop.450 <- round(100/(1+exp(conf.m1[6,]*(log(50000)-log(4000)))),2)
sub.m1.pop.450 <- c(pi.m1.pop.450,ci.m1.pop.450[2],ci.m1.pop.450[1])

pi.m1.claw.10 <- round(100/(1+exp(coef.m1[7]*(1-0))),2)
ci.m1.claw.10 <- round(100/(1+exp(conf.m1[7,]*(1-0))),2)
sub.m1.claw.10 <- c(pi.m1.claw.10,ci.m1.claw.10[2],ci.m1.claw.10[1])

pi.m1.milex.10 <- round(100/(1+exp(coef.m1[3]*(.1-.02))),2)
ci.m1.milex.10 <- round(100/(1+exp(conf.m1[3,]*(.1-.02))),2)
sub.m1.milex.10 <- c(pi.m1.milex.10,ci.m1.milex.10[2],ci.m1.milex.10[1])

pi.m1.milex.20 <- round(100/(1+exp(coef.m1[3]*(.2-.02))),2)
ci.m1.milex.20 <- round(100/(1+exp(conf.m1[3,]*(.2-.02))),2)
sub.m1.milex.20 <- c(pi.m1.milex.20,ci.m1.milex.20[2],ci.m1.milex.20[1])

pi.m1.milex.30 <- round(100/(1+exp(coef.m1[3]*(.3-.02))),2)
ci.m1.milex.30 <- round(100/(1+exp(conf.m1[3,]*(.3-.02))),2)
sub.m1.milex.30 <- c(pi.m1.milex.30,ci.m1.milex.30[2],ci.m1.milex.30[1])


a1.m1 <- rbind(sub.m1.mida.10,sub.m1.mida.20,sub.m1.mida.30)
a2.m1 <- rbind(sub.m1.midb.10,sub.m1.midb.20,sub.m1.midb.30)
a3.m1 <- rbind(sub.m1.milex.10,sub.m1.milex.20,sub.m1.milex.30)
a4.m1 <- rbind(sub.m1.poly.24,sub.m1.poly.26,sub.m1.poly.28)
a5.m1 <- rbind(sub.m1.gdp.24,sub.m1.gdp.26,sub.m1.gdp.212)
a6.m1 <- rbind(sub.m1.pop.410,sub.m1.pop.420,sub.m1.pop.450)
a7.m1 <- sub.m1.claw.10

sub.m1 <- rbind(a1.m1,a2.m1,a3.m1,a4.m1,a5.m1,a6.m1,a7.m1)

## m2
pi.m2.mida.10 <- round(100/(1+exp(coef.m2[1]*(1-0))),2)
ci.m2.mida.10 <- round(100/(1+exp(conf.m2[1,]*(1-0))),2)
sub.m2.mida.10 <- c(pi.m2.mida.10,ci.m2.mida.10[2],ci.m2.mida.10[1])

pi.m2.mida.20 <- round(100/(1+exp(coef.m2[1]*(2-0))),2)
ci.m2.mida.20 <- round(100/(1+exp(conf.m2[1,]*(2-0))),2)
sub.m2.mida.20 <- c(pi.m2.mida.20,ci.m2.mida.20[2],ci.m2.mida.20[1])

pi.m2.mida.30 <- round(100/(1+exp(coef.m2[1]*(3-0))),2)
ci.m2.mida.30 <- round(100/(1+exp(conf.m2[1,]*(3-0))),2)
sub.m2.mida.30 <- c(pi.m2.mida.30,ci.m2.mida.30[2],ci.m2.mida.30[1])

pi.m2.midb.10 <- round(100/(1+exp(coef.m2[2]*(1-0))),2)
ci.m2.midb.10 <- round(100/(1+exp(conf.m2[2,]*(1-0))),2)
sub.m2.midb.10 <- c(pi.m2.midb.10,ci.m2.midb.10[2],ci.m2.midb.10[1])

pi.m2.midb.20 <- round(100/(1+exp(coef.m2[2]*(2-0))),2)
ci.m2.midb.20 <- round(100/(1+exp(conf.m2[2,]*(2-0))),2)
sub.m2.midb.20 <- c(pi.m2.midb.20,ci.m2.midb.20[2],ci.m2.midb.20[1])

pi.m2.midb.30 <- round(100/(1+exp(coef.m2[2]*(3-0))),2)
ci.m2.midb.30 <- round(100/(1+exp(conf.m2[2,]*(3-0))),2)
sub.m2.midb.30 <- c(pi.m2.midb.30,ci.m2.midb.30[2],ci.m2.midb.30[1])

pi.m2.poly.24 <- round(100/(1+exp(coef.m2[4]*(4-2))),2)
ci.m2.poly.24 <- round(100/(1+exp(conf.m2[4,]*(4-2))),2)
sub.m2.poly.24 <- c(pi.m2.poly.24,ci.m2.poly.24[2],ci.m2.poly.24[1])

pi.m2.poly.26 <- round(100/(1+exp(coef.m2[4]*(6-2))),2)
ci.m2.poly.26 <- round(100/(1+exp(conf.m2[4,]*(6-2))),2)
sub.m2.poly.26 <- c(pi.m2.poly.26,ci.m2.poly.26[2],ci.m2.poly.26[1])

pi.m2.poly.28 <- round(100/(1+exp(coef.m2[4]*(8-2))),2)
ci.m2.poly.28 <- round(100/(1+exp(conf.m2[4,]*(8-2))),2)
sub.m2.poly.28 <- c(pi.m2.poly.28,ci.m2.poly.28[2],ci.m2.poly.28[1])

pi.m2.gdp.24 <- round(100/(1+exp(coef.m2[5]*(log(4000)-log(2000)))),2)
ci.m2.gdp.24 <- round(100/(1+exp(conf.m2[5,]*(log(4000)-log(2000)))),2)
sub.m2.gdp.24 <- c(pi.m2.gdp.24,ci.m2.gdp.24[2],ci.m2.gdp.24[1])

pi.m2.gdp.26 <- round(100/(1+exp(coef.m2[5]*(log(8000)-log(2000)))),2)
ci.m2.gdp.26 <- round(100/(1+exp(conf.m2[5,]*(log(8000)-log(2000)))),2)
sub.m2.gdp.26 <- c(pi.m2.gdp.26,ci.m2.gdp.26[2],ci.m2.gdp.26[1])

pi.m2.gdp.212 <- round(100/(1+exp(coef.m2[5]*(log(12000)-log(2000)))),2)
ci.m2.gdp.212 <- round(100/(1+exp(conf.m2[5,]*(log(12000)-log(2000)))),2)
sub.m2.gdp.212 <- c(pi.m2.gdp.212,ci.m2.gdp.212[2],ci.m2.gdp.212[1])

pi.m2.pop.410 <- round(100/(1+exp(coef.m2[6]*(log(10000)-log(4000)))),2)
ci.m2.pop.410 <- round(100/(1+exp(conf.m2[6,]*(log(10000)-log(4000)))),2)
sub.m2.pop.410 <- c(pi.m2.pop.410,ci.m2.pop.410[2],ci.m2.pop.410[1])

pi.m2.pop.420 <- round(100/(1+exp(coef.m2[6]*(log(20000)-log(4000)))),2)
ci.m2.pop.420 <- round(100/(1+exp(conf.m2[6,]*(log(20000)-log(4000)))),2)
sub.m2.pop.420 <- c(pi.m2.pop.420,ci.m2.pop.420[2],ci.m2.pop.420[1])

pi.m2.pop.450 <- round(100/(1+exp(coef.m2[6]*(log(50000)-log(4000)))),2)
ci.m2.pop.450 <- round(100/(1+exp(conf.m2[6,]*(log(50000)-log(4000)))),2)
sub.m2.pop.450 <- c(pi.m2.pop.450,ci.m2.pop.450[2],ci.m2.pop.450[1])

pi.m2.claw.10 <- round(100/(1+exp(coef.m2[7]*(1-0))),2)
ci.m2.claw.10 <- round(100/(1+exp(conf.m2[7,]*(1-0))),2)
sub.m2.claw.10 <- c(pi.m2.claw.10,ci.m2.claw.10[2],ci.m2.claw.10[1])


pi.m2.milex.10 <- round(100/(1+exp(coef.m2[3]*(.1-.02))),2)
ci.m2.milex.10 <- round(100/(1+exp(conf.m2[3,]*(.1-.02))),2)
sub.m2.milex.10 <- c(pi.m2.milex.10,ci.m2.milex.10[2],ci.m2.milex.10[1])

pi.m2.milex.20 <- round(100/(1+exp(coef.m2[3]*(.2-.02))),2)
ci.m2.milex.20 <- round(100/(1+exp(conf.m2[3,]*(.2-.02))),2)
sub.m2.milex.20 <- c(pi.m2.milex.20,ci.m2.milex.20[2],ci.m2.milex.20[1])

pi.m2.milex.30 <- round(100/(1+exp(coef.m2[3]*(.3-.02))),2)
ci.m2.milex.30 <- round(100/(1+exp(conf.m2[3,]*(.3-.02))),2)
sub.m2.milex.30 <- c(pi.m2.milex.30,ci.m2.milex.30[2],ci.m2.milex.30[1])


a1.m2 <- rbind(sub.m2.mida.10,sub.m2.mida.20,sub.m2.mida.30)
a2.m2 <- rbind(sub.m2.midb.10,sub.m2.midb.20,sub.m2.midb.30)
a3.m2 <- rbind(sub.m2.milex.10,sub.m2.milex.20,sub.m2.milex.30)
a4.m2 <- rbind(sub.m2.poly.24,sub.m2.poly.26,sub.m2.poly.28)
a5.m2 <- rbind(sub.m2.gdp.24,sub.m2.gdp.26,sub.m2.gdp.212)
a6.m2 <- rbind(sub.m2.pop.410,sub.m2.pop.420,sub.m2.pop.450)
a7.m2 <- sub.m2.claw.10

sub.m2 <- rbind(a1.m2,a2.m2,a3.m2,a4.m2,a5.m2,a6.m2,a7.m2)


## m3
pi.m3.mida.10 <- round(100/(1+exp(coef.m3[1]*(1-0))),2)
ci.m3.mida.10 <- round(100/(1+exp(conf.m3[1,]*(1-0))),2)
sub.m3.mida.10 <- c(pi.m3.mida.10,ci.m3.mida.10[2],ci.m3.mida.10[1])

pi.m3.mida.20 <- round(100/(1+exp(coef.m3[1]*(2-0))),2)
ci.m3.mida.20 <- round(100/(1+exp(conf.m3[1,]*(2-0))),2)
sub.m3.mida.20 <- c(pi.m3.mida.20,ci.m3.mida.20[2],ci.m3.mida.20[1])

pi.m3.mida.30 <- round(100/(1+exp(coef.m3[1]*(3-0))),2)
ci.m3.mida.30 <- round(100/(1+exp(conf.m3[1,]*(3-0))),2)
sub.m3.mida.30 <- c(pi.m3.mida.30,ci.m3.mida.30[2],ci.m3.mida.30[1])

pi.m3.midb.10 <- round(100/(1+exp(coef.m3[2]*(1-0))),2)
ci.m3.midb.10 <- round(100/(1+exp(conf.m3[2,]*(1-0))),2)
sub.m3.midb.10 <- c(pi.m3.midb.10,ci.m3.midb.10[2],ci.m3.midb.10[1])

pi.m3.midb.20 <- round(100/(1+exp(coef.m3[2]*(2-0))),2)
ci.m3.midb.20 <- round(100/(1+exp(conf.m3[2,]*(2-0))),2)
sub.m3.midb.20 <- c(pi.m3.midb.20,ci.m3.midb.20[2],ci.m3.midb.20[1])

pi.m3.midb.30 <- round(100/(1+exp(coef.m3[2]*(3-0))),2)
ci.m3.midb.30 <- round(100/(1+exp(conf.m3[2,]*(3-0))),2)
sub.m3.midb.30 <- c(pi.m3.midb.30,ci.m3.midb.30[2],ci.m3.midb.30[1])

pi.m3.poly.24 <- round(100/(1+exp(coef.m3[4]*(4-2))),2)
ci.m3.poly.24 <- round(100/(1+exp(conf.m3[4,]*(4-2))),2)
sub.m3.poly.24 <- c(pi.m3.poly.24,ci.m3.poly.24[2],ci.m3.poly.24[1])

pi.m3.poly.26 <- round(100/(1+exp(coef.m3[4]*(6-2))),2)
ci.m3.poly.26 <- round(100/(1+exp(conf.m3[4,]*(6-2))),2)
sub.m3.poly.26 <- c(pi.m3.poly.26,ci.m3.poly.26[2],ci.m3.poly.26[1])

pi.m3.poly.28 <- round(100/(1+exp(coef.m3[4]*(8-2))),2)
ci.m3.poly.28 <- round(100/(1+exp(conf.m3[4,]*(8-2))),2)
sub.m3.poly.28 <- c(pi.m3.poly.28,ci.m3.poly.28[2],ci.m3.poly.28[1])

pi.m3.gdp.24 <- round(100/(1+exp(coef.m3[5]*(log(4000)-log(2000)))),2)
ci.m3.gdp.24 <- round(100/(1+exp(conf.m3[5,]*(log(4000)-log(2000)))),2)
sub.m3.gdp.24 <- c(pi.m3.gdp.24,ci.m3.gdp.24[2],ci.m3.gdp.24[1])

pi.m3.gdp.26 <- round(100/(1+exp(coef.m3[5]*(log(8000)-log(2000)))),2)
ci.m3.gdp.26 <- round(100/(1+exp(conf.m3[5,]*(log(8000)-log(2000)))),2)
sub.m3.gdp.26 <- c(pi.m3.gdp.26,ci.m3.gdp.26[2],ci.m3.gdp.26[1])

pi.m3.gdp.212 <- round(100/(1+exp(coef.m3[5]*(log(12000)-log(2000)))),2)
ci.m3.gdp.212 <- round(100/(1+exp(conf.m3[5,]*(log(12000)-log(2000)))),2)
sub.m3.gdp.212 <- c(pi.m3.gdp.212,ci.m3.gdp.212[2],ci.m3.gdp.212[1])

pi.m3.pop.410 <- round(100/(1+exp(coef.m3[6]*(log(10000)-log(4000)))),2)
ci.m3.pop.410 <- round(100/(1+exp(conf.m3[6,]*(log(10000)-log(4000)))),2)
sub.m3.pop.410 <- c(pi.m3.pop.410,ci.m3.pop.410[2],ci.m3.pop.410[1])

pi.m3.pop.420 <- round(100/(1+exp(coef.m3[6]*(log(20000)-log(4000)))),2)
ci.m3.pop.420 <- round(100/(1+exp(conf.m3[6,]*(log(20000)-log(4000)))),2)
sub.m3.pop.420 <- c(pi.m3.pop.420,ci.m3.pop.420[2],ci.m3.pop.420[1])

pi.m3.pop.450 <- round(100/(1+exp(coef.m3[6]*(log(50000)-log(4000)))),2)
ci.m3.pop.450 <- round(100/(1+exp(conf.m3[6,]*(log(50000)-log(4000)))),2)
sub.m3.pop.450 <- c(pi.m3.pop.450,ci.m3.pop.450[2],ci.m3.pop.450[1])

pi.m3.claw.10 <- round(100/(1+exp(coef.m3[7]*(1-0))),2)
ci.m3.claw.10 <- round(100/(1+exp(conf.m3[7,]*(1-0))),2)
sub.m3.claw.10 <- c(pi.m3.claw.10,ci.m3.claw.10[2],ci.m3.claw.10[1])

pi.m3.milex.10 <- round(100/(1+exp(coef.m3[3]*(.1-.02))),2)
ci.m3.milex.10 <- round(100/(1+exp(conf.m3[3,]*(.1-.02))),2)
sub.m3.milex.10 <- c(pi.m3.milex.10,ci.m3.milex.10[2],ci.m3.milex.10[1])

pi.m3.milex.20 <- round(100/(1+exp(coef.m3[3]*(.2-.02))),2)
ci.m3.milex.20 <- round(100/(1+exp(conf.m3[3,]*(.2-.02))),2)
sub.m3.milex.20 <- c(pi.m3.milex.20,ci.m3.milex.20[2],ci.m3.milex.20[1])

pi.m3.milex.30 <- round(100/(1+exp(coef.m3[3]*(.3-.02))),2)
ci.m3.milex.30 <- round(100/(1+exp(conf.m3[3,]*(.3-.02))),2)
sub.m3.milex.30 <- c(pi.m3.milex.30,ci.m3.milex.30[2],ci.m3.milex.30[1])


a1.m3 <- rbind(sub.m3.mida.10,sub.m3.mida.20,sub.m3.mida.30)
a2.m3 <- rbind(sub.m3.midb.10,sub.m3.midb.20,sub.m3.midb.30)
a3.m3 <- rbind(sub.m3.milex.10,sub.m3.milex.20,sub.m3.milex.30)
a4.m3 <- rbind(sub.m3.poly.24,sub.m3.poly.26,sub.m3.poly.28)
a5.m3 <- rbind(sub.m3.gdp.24,sub.m3.gdp.26,sub.m3.gdp.212)
a6.m3 <- rbind(sub.m3.pop.410,sub.m3.pop.420,sub.m3.pop.450)
a7.m3 <- sub.m3.claw.10

sub.m3 <- rbind(a1.m3,a2.m3,a3.m3,a4.m3,a5.m3,a6.m3,a7.m3)

## m4
pi.m4.mida.10 <- round(100/(1+exp(coef.m4[1]*(1-0))),2)
ci.m4.mida.10 <- round(100/(1+exp(conf.m4[1,]*(1-0))),2)
sub.m4.mida.10 <- c(pi.m4.mida.10,ci.m4.mida.10[2],ci.m4.mida.10[1])

pi.m4.mida.20 <- round(100/(1+exp(coef.m4[1]*(2-0))),2)
ci.m4.mida.20 <- round(100/(1+exp(conf.m4[1,]*(2-0))),2)
sub.m4.mida.20 <- c(pi.m4.mida.20,ci.m4.mida.20[2],ci.m4.mida.20[1])

pi.m4.mida.30 <- round(100/(1+exp(coef.m4[1]*(3-0))),2)
ci.m4.mida.30 <- round(100/(1+exp(conf.m4[1,]*(3-0))),2)
sub.m4.mida.30 <- c(pi.m4.mida.30,ci.m4.mida.30[2],ci.m4.mida.30[1])

pi.m4.midb.10 <- round(100/(1+exp(coef.m4[2]*(1-0))),2)
ci.m4.midb.10 <- round(100/(1+exp(conf.m4[2,]*(1-0))),2)
sub.m4.midb.10 <- c(pi.m4.midb.10,ci.m4.midb.10[2],ci.m4.midb.10[1])

pi.m4.midb.20 <- round(100/(1+exp(coef.m4[2]*(2-0))),2)
ci.m4.midb.20 <- round(100/(1+exp(conf.m4[2,]*(2-0))),2)
sub.m4.midb.20 <- c(pi.m4.midb.20,ci.m4.midb.20[2],ci.m4.midb.20[1])

pi.m4.midb.30 <- round(100/(1+exp(coef.m4[2]*(3-0))),2)
ci.m4.midb.30 <- round(100/(1+exp(conf.m4[2,]*(3-0))),2)
sub.m4.midb.30 <- c(pi.m4.midb.30,ci.m4.midb.30[2],ci.m4.midb.30[1])

pi.m4.poly.24 <- round(100/(1+exp(coef.m4[4]*(4-2))),2)
ci.m4.poly.24 <- round(100/(1+exp(conf.m4[4,]*(4-2))),2)
sub.m4.poly.24 <- c(pi.m4.poly.24,ci.m4.poly.24[2],ci.m4.poly.24[1])

pi.m4.poly.26 <- round(100/(1+exp(coef.m4[4]*(6-2))),2)
ci.m4.poly.26 <- round(100/(1+exp(conf.m4[4,]*(6-2))),2)
sub.m4.poly.26 <- c(pi.m4.poly.26,ci.m4.poly.26[2],ci.m4.poly.26[1])

pi.m4.poly.28 <- round(100/(1+exp(coef.m4[4]*(8-2))),2)
ci.m4.poly.28 <- round(100/(1+exp(conf.m4[4,]*(8-2))),2)
sub.m4.poly.28 <- c(pi.m4.poly.28,ci.m4.poly.28[2],ci.m4.poly.28[1])

pi.m4.gdp.24 <- round(100/(1+exp(coef.m4[5]*(log(4000)-log(2000)))),2)
ci.m4.gdp.24 <- round(100/(1+exp(conf.m4[5,]*(log(4000)-log(2000)))),2)
sub.m4.gdp.24 <- c(pi.m4.gdp.24,ci.m4.gdp.24[2],ci.m4.gdp.24[1])

pi.m4.gdp.26 <- round(100/(1+exp(coef.m4[5]*(log(8000)-log(2000)))),2)
ci.m4.gdp.26 <- round(100/(1+exp(conf.m4[5,]*(log(8000)-log(2000)))),2)
sub.m4.gdp.26 <- c(pi.m4.gdp.26,ci.m4.gdp.26[2],ci.m4.gdp.26[1])

pi.m4.gdp.212 <- round(100/(1+exp(coef.m4[5]*(log(12000)-log(2000)))),2)
ci.m4.gdp.212 <- round(100/(1+exp(conf.m4[5,]*(log(12000)-log(2000)))),2)
sub.m4.gdp.212 <- c(pi.m4.gdp.212,ci.m4.gdp.212[2],ci.m4.gdp.212[1])

pi.m4.pop.410 <- round(100/(1+exp(coef.m4[6]*(log(10000)-log(4000)))),2)
ci.m4.pop.410 <- round(100/(1+exp(conf.m4[6,]*(log(10000)-log(4000)))),2)
sub.m4.pop.410 <- c(pi.m4.pop.410,ci.m4.pop.410[2],ci.m4.pop.410[1])

pi.m4.pop.420 <- round(100/(1+exp(coef.m4[6]*(log(20000)-log(4000)))),2)
ci.m4.pop.420 <- round(100/(1+exp(conf.m4[6,]*(log(20000)-log(4000)))),2)
sub.m4.pop.420 <- c(pi.m4.pop.420,ci.m4.pop.420[2],ci.m4.pop.420[1])

pi.m4.pop.450 <- round(100/(1+exp(coef.m4[6]*(log(50000)-log(4000)))),2)
ci.m4.pop.450 <- round(100/(1+exp(conf.m4[6,]*(log(50000)-log(4000)))),2)
sub.m4.pop.450 <- c(pi.m4.pop.450,ci.m4.pop.450[2],ci.m4.pop.450[1])

pi.m4.claw.10 <- round(100/(1+exp(coef.m4[7]*(1-0))),2)
ci.m4.claw.10 <- round(100/(1+exp(conf.m4[7,]*(1-0))),2)
sub.m4.claw.10 <- c(pi.m4.claw.10,ci.m4.claw.10[2],ci.m4.claw.10[1])

pi.m4.milex.10 <- round(100/(1+exp(coef.m4[3]*(.1-.02))),2)
ci.m4.milex.10 <- round(100/(1+exp(conf.m4[3,]*(.1-.02))),2)
sub.m4.milex.10 <- c(pi.m4.milex.10,ci.m4.milex.10[2],ci.m4.milex.10[1])

pi.m4.milex.20 <- round(100/(1+exp(coef.m4[3]*(.2-.02))),2)
ci.m4.milex.20 <- round(100/(1+exp(conf.m4[3,]*(.2-.02))),2)
sub.m4.milex.20 <- c(pi.m4.milex.20,ci.m4.milex.20[2],ci.m4.milex.20[1])

pi.m4.milex.30 <- round(100/(1+exp(coef.m4[3]*(.3-.02))),2)
ci.m4.milex.30 <- round(100/(1+exp(conf.m4[3,]*(.3-.02))),2)
sub.m4.milex.30 <- c(pi.m4.milex.30,ci.m4.milex.30[2],ci.m4.milex.30[1])


a1.m4 <- rbind(sub.m4.mida.10,sub.m4.mida.20,sub.m4.mida.30)
a2.m4 <- rbind(sub.m4.midb.10,sub.m4.midb.20,sub.m4.midb.30)
a3.m4 <- rbind(sub.m4.milex.10,sub.m4.milex.20,sub.m4.milex.30)
a4.m4 <- rbind(sub.m4.poly.24,sub.m4.poly.26,sub.m4.poly.28)
a5.m4 <- rbind(sub.m4.gdp.24,sub.m4.gdp.26,sub.m4.gdp.212)
a6.m4 <- rbind(sub.m4.pop.410,sub.m4.pop.420,sub.m4.pop.450)
a7.m4 <- sub.m4.claw.10

sub.m4 <- rbind(a1.m4,a2.m4,a3.m4,a4.m4,a5.m4,a6.m4,a7.m4)


## m5
pi.m5.mida.10 <- round(100/(1+exp(coef.m5[1]*(1-0))),2)
ci.m5.mida.10 <- round(100/(1+exp(conf.m5[1,]*(1-0))),2)
sub.m5.mida.10 <- c(pi.m5.mida.10,ci.m5.mida.10[2],ci.m5.mida.10[1])

pi.m5.mida.20 <- round(100/(1+exp(coef.m5[1]*(2-0))),2)
ci.m5.mida.20 <- round(100/(1+exp(conf.m5[1,]*(2-0))),2)
sub.m5.mida.20 <- c(pi.m5.mida.20,ci.m5.mida.20[2],ci.m5.mida.20[1])

pi.m5.mida.30 <- round(100/(1+exp(coef.m5[1]*(3-0))),2)
ci.m5.mida.30 <- round(100/(1+exp(conf.m5[1,]*(3-0))),2)
sub.m5.mida.30 <- c(pi.m5.mida.30,ci.m5.mida.30[2],ci.m5.mida.30[1])

pi.m5.midb.10 <- round(100/(1+exp(coef.m5[2]*(1-0))),2)
ci.m5.midb.10 <- round(100/(1+exp(conf.m5[2,]*(1-0))),2)
sub.m5.midb.10 <- c(pi.m5.midb.10,ci.m5.midb.10[2],ci.m5.midb.10[1])

pi.m5.midb.20 <- round(100/(1+exp(coef.m5[2]*(2-0))),2)
ci.m5.midb.20 <- round(100/(1+exp(conf.m5[2,]*(2-0))),2)
sub.m5.midb.20 <- c(pi.m5.midb.20,ci.m5.midb.20[2],ci.m5.midb.20[1])

pi.m5.midb.30 <- round(100/(1+exp(coef.m5[2]*(3-0))),2)
ci.m5.midb.30 <- round(100/(1+exp(conf.m5[2,]*(3-0))),2)
sub.m5.midb.30 <- c(pi.m5.midb.30,ci.m5.midb.30[2],ci.m5.midb.30[1])

pi.m5.poly.24 <- round(100/(1+exp(coef.m5[4]*(4-2))),2)
ci.m5.poly.24 <- round(100/(1+exp(conf.m5[4,]*(4-2))),2)
sub.m5.poly.24 <- c(pi.m5.poly.24,ci.m5.poly.24[2],ci.m5.poly.24[1])

pi.m5.poly.26 <- round(100/(1+exp(coef.m5[4]*(6-2))),2)
ci.m5.poly.26 <- round(100/(1+exp(conf.m5[4,]*(6-2))),2)
sub.m5.poly.26 <- c(pi.m5.poly.26,ci.m5.poly.26[2],ci.m5.poly.26[1])

pi.m5.poly.28 <- round(100/(1+exp(coef.m5[4]*(8-2))),2)
ci.m5.poly.28 <- round(100/(1+exp(conf.m5[4,]*(8-2))),2)
sub.m5.poly.28 <- c(pi.m5.poly.28,ci.m5.poly.28[2],ci.m5.poly.28[1])

pi.m5.gdp.24 <- round(100/(1+exp(coef.m5[5]*(log(4000)-log(2000)))),2)
ci.m5.gdp.24 <- round(100/(1+exp(conf.m5[5,]*(log(4000)-log(2000)))),2)
sub.m5.gdp.24 <- c(pi.m5.gdp.24,ci.m5.gdp.24[2],ci.m5.gdp.24[1])

pi.m5.gdp.26 <- round(100/(1+exp(coef.m5[5]*(log(8000)-log(2000)))),2)
ci.m5.gdp.26 <- round(100/(1+exp(conf.m5[5,]*(log(8000)-log(2000)))),2)
sub.m5.gdp.26 <- c(pi.m5.gdp.26,ci.m5.gdp.26[2],ci.m5.gdp.26[1])

pi.m5.gdp.212 <- round(100/(1+exp(coef.m5[5]*(log(12000)-log(2000)))),2)
ci.m5.gdp.212 <- round(100/(1+exp(conf.m5[5,]*(log(12000)-log(2000)))),2)
sub.m5.gdp.212 <- c(pi.m5.gdp.212,ci.m5.gdp.212[2],ci.m5.gdp.212[1])

pi.m5.pop.410 <- round(100/(1+exp(coef.m5[6]*(log(10000)-log(4000)))),2)
ci.m5.pop.410 <- round(100/(1+exp(conf.m5[6,]*(log(10000)-log(4000)))),2)
sub.m5.pop.410 <- c(pi.m5.pop.410,ci.m5.pop.410[2],ci.m5.pop.410[1])

pi.m5.pop.420 <- round(100/(1+exp(coef.m5[6]*(log(20000)-log(4000)))),2)
ci.m5.pop.420 <- round(100/(1+exp(conf.m5[6,]*(log(20000)-log(4000)))),2)
sub.m5.pop.420 <- c(pi.m5.pop.420,ci.m5.pop.420[2],ci.m5.pop.420[1])

pi.m5.pop.450 <- round(100/(1+exp(coef.m5[6]*(log(50000)-log(4000)))),2)
ci.m5.pop.450 <- round(100/(1+exp(conf.m5[6,]*(log(50000)-log(4000)))),2)
sub.m5.pop.450 <- c(pi.m5.pop.450,ci.m5.pop.450[2],ci.m5.pop.450[1])

pi.m5.claw.10 <- round(100/(1+exp(coef.m5[7]*(1-0))),2)
ci.m5.claw.10 <- round(100/(1+exp(conf.m5[7,]*(1-0))),2)
sub.m5.claw.10 <- c(pi.m5.claw.10,ci.m5.claw.10[2],ci.m5.claw.10[1])

pi.m5.milex.10 <- round(100/(1+exp(coef.m5[3]*(.1-.02))),2)
ci.m5.milex.10 <- round(100/(1+exp(conf.m5[3,]*(.1-.02))),2)
sub.m5.milex.10 <- c(pi.m5.milex.10,ci.m5.milex.10[2],ci.m5.milex.10[1])

pi.m5.milex.20 <- round(100/(1+exp(coef.m5[3]*(.2-.02))),2)
ci.m5.milex.20 <- round(100/(1+exp(conf.m5[3,]*(.2-.02))),2)
sub.m5.milex.20 <- c(pi.m5.milex.20,ci.m5.milex.20[2],ci.m5.milex.20[1])

pi.m5.milex.30 <- round(100/(1+exp(coef.m5[3]*(.3-.02))),2)
ci.m5.milex.30 <- round(100/(1+exp(conf.m5[3,]*(.3-.02))),2)
sub.m5.milex.30 <- c(pi.m5.milex.30,ci.m5.milex.30[2],ci.m5.milex.30[1])


a1.m5 <- rbind(sub.m5.mida.10,sub.m5.mida.20,sub.m5.mida.30)
a2.m5 <- rbind(sub.m5.midb.10,sub.m5.midb.20,sub.m5.midb.30)
a3.m5 <- rbind(sub.m5.milex.10,sub.m5.milex.20,sub.m5.milex.30)
a4.m5 <- rbind(sub.m5.poly.24,sub.m5.poly.26,sub.m5.poly.28)
a5.m5 <- rbind(sub.m5.gdp.24,sub.m5.gdp.26,sub.m5.gdp.212)
a6.m5 <- rbind(sub.m5.pop.410,sub.m5.pop.420,sub.m5.pop.450)
a7.m5 <- sub.m5.claw.10

sub.m5 <- rbind(a1.m5,a2.m5,a3.m5,a4.m5,a5.m5,a6.m5,a7.m5)


## m6
pi.m6.mida.10 <- round(100/(1+exp(coef.m6[1]*(1-0))),2)
ci.m6.mida.10 <- round(100/(1+exp(conf.m6[1,]*(1-0))),2)
sub.m6.mida.10 <- c(pi.m6.mida.10,ci.m6.mida.10[2],ci.m6.mida.10[1])

pi.m6.mida.20 <- round(100/(1+exp(coef.m6[1]*(2-0))),2)
ci.m6.mida.20 <- round(100/(1+exp(conf.m6[1,]*(2-0))),2)
sub.m6.mida.20 <- c(pi.m6.mida.20,ci.m6.mida.20[2],ci.m6.mida.20[1])

pi.m6.mida.30 <- round(100/(1+exp(coef.m6[1]*(3-0))),2)
ci.m6.mida.30 <- round(100/(1+exp(conf.m6[1,]*(3-0))),2)
sub.m6.mida.30 <- c(pi.m6.mida.30,ci.m6.mida.30[2],ci.m6.mida.30[1])

pi.m6.midb.10 <- round(100/(1+exp(coef.m6[2]*(1-0))),2)
ci.m6.midb.10 <- round(100/(1+exp(conf.m6[2,]*(1-0))),2)
sub.m6.midb.10 <- c(pi.m6.midb.10,ci.m6.midb.10[2],ci.m6.midb.10[1])

pi.m6.midb.20 <- round(100/(1+exp(coef.m6[2]*(2-0))),2)
ci.m6.midb.20 <- round(100/(1+exp(conf.m6[2,]*(2-0))),2)
sub.m6.midb.20 <- c(pi.m6.midb.20,ci.m6.midb.20[2],ci.m6.midb.20[1])

pi.m6.midb.30 <- round(100/(1+exp(coef.m6[2]*(3-0))),2)
ci.m6.midb.30 <- round(100/(1+exp(conf.m6[2,]*(3-0))),2)
sub.m6.midb.30 <- c(pi.m6.midb.30,ci.m6.midb.30[2],ci.m6.midb.30[1])

pi.m6.poly.24 <- round(100/(1+exp(coef.m6[4]*(4-2))),2)
ci.m6.poly.24 <- round(100/(1+exp(conf.m6[4,]*(4-2))),2)
sub.m6.poly.24 <- c(pi.m6.poly.24,ci.m6.poly.24[2],ci.m6.poly.24[1])

pi.m6.poly.26 <- round(100/(1+exp(coef.m6[4]*(6-2))),2)
ci.m6.poly.26 <- round(100/(1+exp(conf.m6[4,]*(6-2))),2)
sub.m6.poly.26 <- c(pi.m6.poly.26,ci.m6.poly.26[2],ci.m6.poly.26[1])

pi.m6.poly.28 <- round(100/(1+exp(coef.m6[4]*(8-2))),2)
ci.m6.poly.28 <- round(100/(1+exp(conf.m6[4,]*(8-2))),2)
sub.m6.poly.28 <- c(pi.m6.poly.28,ci.m6.poly.28[2],ci.m6.poly.28[1])

pi.m6.gdp.24 <- round(100/(1+exp(coef.m6[5]*(log(4000)-log(2000)))),2)
ci.m6.gdp.24 <- round(100/(1+exp(conf.m6[5,]*(log(4000)-log(2000)))),2)
sub.m6.gdp.24 <- c(pi.m6.gdp.24,ci.m6.gdp.24[2],ci.m6.gdp.24[1])

pi.m6.gdp.26 <- round(100/(1+exp(coef.m6[5]*(log(8000)-log(2000)))),2)
ci.m6.gdp.26 <- round(100/(1+exp(conf.m6[5,]*(log(8000)-log(2000)))),2)
sub.m6.gdp.26 <- c(pi.m6.gdp.26,ci.m6.gdp.26[2],ci.m6.gdp.26[1])

pi.m6.gdp.212 <- round(100/(1+exp(coef.m6[5]*(log(12000)-log(2000)))),2)
ci.m6.gdp.212 <- round(100/(1+exp(conf.m6[5,]*(log(12000)-log(2000)))),2)
sub.m6.gdp.212 <- c(pi.m6.gdp.212,ci.m6.gdp.212[2],ci.m6.gdp.212[1])

pi.m6.pop.410 <- round(100/(1+exp(coef.m6[6]*(log(10000)-log(4000)))),2)
ci.m6.pop.410 <- round(100/(1+exp(conf.m6[6,]*(log(10000)-log(4000)))),2)
sub.m6.pop.410 <- c(pi.m6.pop.410,ci.m6.pop.410[2],ci.m6.pop.410[1])

pi.m6.pop.420 <- round(100/(1+exp(coef.m6[6]*(log(20000)-log(4000)))),2)
ci.m6.pop.420 <- round(100/(1+exp(conf.m6[6,]*(log(20000)-log(4000)))),2)
sub.m6.pop.420 <- c(pi.m6.pop.420,ci.m6.pop.420[2],ci.m6.pop.420[1])

pi.m6.pop.450 <- round(100/(1+exp(coef.m6[6]*(log(50000)-log(4000)))),2)
ci.m6.pop.450 <- round(100/(1+exp(conf.m6[6,]*(log(50000)-log(4000)))),2)
sub.m6.pop.450 <- c(pi.m6.pop.450,ci.m6.pop.450[2],ci.m6.pop.450[1])

pi.m6.claw.10 <- round(100/(1+exp(coef.m6[7]*(1-0))),2)
ci.m6.claw.10 <- round(100/(1+exp(conf.m6[7,]*(1-0))),2)
sub.m6.claw.10 <- c(pi.m6.claw.10,ci.m6.claw.10[2],ci.m6.claw.10[1])

pi.m6.milex.10 <- round(100/(1+exp(coef.m6[3]*(.1-.02))),2)
ci.m6.milex.10 <- round(100/(1+exp(conf.m6[3,]*(.1-.02))),2)
sub.m6.milex.10 <- c(pi.m6.milex.10,ci.m6.milex.10[2],ci.m6.milex.10[1])

pi.m6.milex.20 <- round(100/(1+exp(coef.m6[3]*(.2-.02))),2)
ci.m6.milex.20 <- round(100/(1+exp(conf.m6[3,]*(.2-.02))),2)
sub.m6.milex.20 <- c(pi.m6.milex.20,ci.m6.milex.20[2],ci.m6.milex.20[1])

pi.m6.milex.30 <- round(100/(1+exp(coef.m6[3]*(.3-.02))),2)
ci.m6.milex.30 <- round(100/(1+exp(conf.m6[3,]*(.3-.02))),2)
sub.m6.milex.30 <- c(pi.m6.milex.30,ci.m6.milex.30[2],ci.m6.milex.30[1])

a1.m6 <- rbind(sub.m6.mida.10,sub.m6.mida.20,sub.m6.mida.30)
a2.m6 <- rbind(sub.m6.midb.10,sub.m6.midb.20,sub.m6.midb.30)
a3.m6 <- rbind(sub.m6.milex.10,sub.m6.milex.20,sub.m6.milex.30)
a4.m6 <- rbind(sub.m6.poly.24,sub.m6.poly.26,sub.m6.poly.28)
a5.m6 <- rbind(sub.m6.gdp.24,sub.m6.gdp.26,sub.m6.gdp.212)
a6.m6 <- rbind(sub.m6.pop.410,sub.m6.pop.420,sub.m6.pop.450)
a7.m6 <- sub.m6.claw.10

sub.m6 <- rbind(a1.m6,a2.m6,a3.m6,a4.m6,a5.m6,a6.m6,a7.m6)

sub.fig1.a <- cbind(sub.m1,sub.m3,sub.m5)
sub.fig1.b <- cbind(sub.m2,sub.m4,sub.m6)

sub.fig1 <- data.frame(rbind(sub.fig1.a,sub.fig1.b))

sub.fig1$contrast <- NULL
sub.fig1$contrast[1] <- "1 vs. 0"
sub.fig1$contrast[2] <- "2 vs. 0"
sub.fig1$contrast[3] <- "3 vs. 0"
sub.fig1$contrast[4] <- "1 vs. 0"
sub.fig1$contrast[5] <- "2 vs. 0"
sub.fig1$contrast[6] <- "3 vs. 0"
sub.fig1$contrast[7] <- "100M vs. 20M"
sub.fig1$contrast[8] <- "200M vs. 20M"
sub.fig1$contrast[9] <- "300M vs. 20M"
sub.fig1$contrast[10] <- "4 vs. 2"
sub.fig1$contrast[11] <- "6 vs. 2"
sub.fig1$contrast[12] <- "8 vs. 2"
sub.fig1$contrast[13] <- "4K vs. 2K"
sub.fig1$contrast[14] <- "8K vs. 2K"
sub.fig1$contrast[15] <- "12K vs. 2K"
sub.fig1$contrast[16] <- "10M vs. 4M"
sub.fig1$contrast[17] <- "20M vs. 4M"
sub.fig1$contrast[18] <- "50M vs. 4M"
sub.fig1$contrast[19] <- "Yes vs. No"
sub.fig1$contrast[20] <- "1 vs. 0"
sub.fig1$contrast[21] <- "2 vs. 0"
sub.fig1$contrast[22] <- "3 vs. 0"
sub.fig1$contrast[23] <- "1 vs. 0"
sub.fig1$contrast[24] <- "2 vs. 0"
sub.fig1$contrast[25] <- "3 vs. 0"
sub.fig1$contrast[26] <- "100M vs. 20M"
sub.fig1$contrast[27] <- "200M vs. 20M"
sub.fig1$contrast[28] <- "300M vs. 20M"
sub.fig1$contrast[29] <- "4 vs. 2"
sub.fig1$contrast[30] <- "6 vs. 2"
sub.fig1$contrast[31] <- "8 vs. 2"
sub.fig1$contrast[32] <- "4K vs. 2K"
sub.fig1$contrast[33] <- "8K vs. 2K"
sub.fig1$contrast[34] <- "12K vs. 2K"
sub.fig1$contrast[35] <- "10M vs. 4M"
sub.fig1$contrast[36] <- "20M vs. 4M"
sub.fig1$contrast[37] <- "50M vs. 4M"
sub.fig1$contrast[38] <- "Yes vs. No"

sub.fig1$contrast1 <- NULL
sub.fig1$contrast1[1] <- "MID side A: 1 vs. 0"
sub.fig1$contrast1[2] <- "2 vs. 0"
sub.fig1$contrast1[3] <- "3 vs. 0"
sub.fig1$contrast1[4] <- "MID side B: 1 vs. 0"
sub.fig1$contrast1[5] <- "2  vs. 0"
sub.fig1$contrast1[6] <- "3  vs. 0"
sub.fig1$contrast1[7] <- "Milex cap: 100M vs. 20M"
sub.fig1$contrast1[8] <- "200M vs. 20M"
sub.fig1$contrast1[9] <- "300M vs. 20M"
sub.fig1$contrast1[10] <- "Polyarchy: 4 vs. 2"
sub.fig1$contrast1[11] <- "6 vs. 2"
sub.fig1$contrast1[12] <- "8 vs. 2"
sub.fig1$contrast1[13] <- "GDP cap (log): 4K vs. 2K"
sub.fig1$contrast1[14] <- "8K vs. 2K"
sub.fig1$contrast1[15] <- "12K vs. 2K"
sub.fig1$contrast1[16] <- "Population: 10M vs. 4M"
sub.fig1$contrast1[17] <- "20M vs. 4M"
sub.fig1$contrast1[18] <- "50M vs. 4M"
sub.fig1$contrast1[19] <- "Common law: Yes vs. No"
sub.fig1$contrast1[20] <- "MID side A: 1 vs. 0"       
sub.fig1$contrast1[21] <- "2 vs. 0"       
sub.fig1$contrast1[22] <- "3 vs. 0"       
sub.fig1$contrast1[23] <- "MID side B: 1 vs. 0"       
sub.fig1$contrast1[24] <- "2  vs. 0"       
sub.fig1$contrast1[25] <- "3  vs. 0"
sub.fig1$contrast1[26] <-  "Milex cap: 100M vs. 20M"
sub.fig1$contrast1[27] <-  "200M vs. 20M"           
sub.fig1$contrast1[28] <-  "300M vs. 20M"           
sub.fig1$contrast1[29] <- "Polyarchy: 4 vs. 2"      
sub.fig1$contrast1[30] <- "6 vs. 2"      
sub.fig1$contrast1[31] <- "8 vs. 2"      
sub.fig1$contrast1[32] <-  "GDP cap (log): 4K vs. 2K" 
sub.fig1$contrast1[33] <-  "8K vs. 2K" 
sub.fig1$contrast1[34] <-  "12K vs. 2K"
sub.fig1$contrast1[35] <-  "Population: 10M vs. 4M"   
sub.fig1$contrast1[36] <-  "20M vs. 4M"   
sub.fig1$contrast1[37] <-  "50M vs. 4M"   
sub.fig1$contrast1[38] <-  "Common law: Yes vs. No"   


sub.fig1$Variables <- NULL
sub.fig1$Variables[1] <- "MID side A"
sub.fig1$Variables[2] <- "MID side A"
sub.fig1$Variables[3] <- "MID side A"
sub.fig1$Variables[4] <- "MID side B"
sub.fig1$Variables[5] <- "MID side B"
sub.fig1$Variables[6] <- "MID side B"
sub.fig1$Variables[7] <- "Milex cap"
sub.fig1$Variables[8] <- "Milex cap"
sub.fig1$Variables[9] <- "Milex cap"
sub.fig1$Variables[10] <- "Polyarchy"
sub.fig1$Variables[11] <- "Polyarchy"
sub.fig1$Variables[12] <- "Polyarchy"
sub.fig1$Variables[13] <- "GDP cap (log)"
sub.fig1$Variables[14] <- "GDP cap (log)"
sub.fig1$Variables[15] <- "GDP cap (log)"
sub.fig1$Variables[16] <- "Population (log)"
sub.fig1$Variables[17] <- "Population (log)"
sub.fig1$Variables[18] <- "Population (log)"
sub.fig1$Variables[19] <- "Common law"
sub.fig1$Variables[20] <- "MID side A"       
sub.fig1$Variables[21] <- "MID side A"       
sub.fig1$Variables[22] <- "MID side A"       
sub.fig1$Variables[23] <- "MID side B"       
sub.fig1$Variables[24] <- "MID side B"       
sub.fig1$Variables[25] <- "MID side B"
sub.fig1$Variables[26] <- "Milex cap"
sub.fig1$Variables[27] <- "Milex cap"
sub.fig1$Variables[28] <- "Milex cap"
sub.fig1$Variables[29] <- "Polyarchy"        
sub.fig1$Variables[30] <- "Polyarchy"        
sub.fig1$Variables[31] <- "Polyarchy"        
sub.fig1$Variables[32] <- "GDP cap (log)"   
sub.fig1$Variables[33] <- "GDP cap (log)"   
sub.fig1$Variables[34] <- "GDP cap (log)"   
sub.fig1$Variables[35] <- "Population (log)"
sub.fig1$Variables[36] <- "Population (log)"
sub.fig1$Variables[37] <- "Population (log)"
sub.fig1$Variables[38] <- "Common law"      


a.fig1 <- sub.fig1[1:19,c(1:3,10:12)]
b.fig1 <- sub.fig1[1:19,c(4:6,10:12)]
c.fig1 <- sub.fig1[1:19,c(7:9,10:12)]
d.fig1 <- sub.fig1[20:38,c(1:3,10:12)]
e.fig1 <- sub.fig1[20:38,c(4:6,10:12)]
f.fig1 <- sub.fig1[20:38,c(7:9,10:12)]

colnames(a.fig1) <- c("PI","X97.5","X2.5","Contrast","Contrast1","Variables")
colnames(b.fig1) <- c("PI","X97.5","X2.5","Contrast","Contrast1","Variables")
colnames(c.fig1) <- c("PI","X97.5","X2.5","Contrast","Contrast1","Variables")
colnames(d.fig1) <- c("PI","X97.5","X2.5","Contrast","Contrast1","Variables")
colnames(e.fig1) <- c("PI","X97.5","X2.5","Contrast","Contrast1","Variables")
colnames(f.fig1) <- c("PI","X97.5","X2.5","Contrast","Contrast1","Variables")

a.fig1$type <- NULL
a.fig1$type <- "Overall,\nfrom 1991"
b.fig1$type <- NULL
b.fig1$type <- "Force,\nfrom 1991"
c.fig1$type <- NULL
c.fig1$type <- "No force,\nfrom 1991"
d.fig1$type <- NULL
d.fig1$type <- "Overall,\nfrom 1945"
e.fig1$type <- NULL
e.fig1$type <- "Force,\nfrom 1945"
f.fig1$type <- NULL
f.fig1$type <- "No force,\nfrom 1945"

a.fig1$type1 <- NULL
a.fig1$type1 <- "Overall"
b.fig1$type1 <- NULL
b.fig1$type1 <- "Force"
c.fig1$type1 <- NULL
c.fig1$type1 <- "No force"
d.fig1$type1 <- NULL
d.fig1$type1 <- "Overall"
e.fig1$type1 <- NULL
e.fig1$type1 <- "Force"
f.fig1$type1 <- NULL
f.fig1$type1 <- "No force"

a.fig1$type2 <- NULL
a.fig1$type2 <- "from 1991"
b.fig1$type2 <- NULL
b.fig1$type2 <- "from 1991"
c.fig1$type2 <- NULL
c.fig1$type2 <- "from 1991"
d.fig1$type2 <- NULL
d.fig1$type2 <- "from 1945"
e.fig1$type2 <- NULL
e.fig1$type2 <- "from 1945"
f.fig1$type2 <- NULL
f.fig1$type2 <- "from 1945"

sub.bigfig <- rbind(a.fig1,
                    b.fig1,
                    c.fig1,
                    d.fig1,
                    e.fig1,
                    f.fig1)

sub.bigfig$type.factor <- factor(sub.bigfig$type,
                                 levels=c("No force,\nfrom 1945",
                                          "Force,\nfrom 1945",
                                          "Overall,\nfrom 1945",
                                          "No force,\nfrom 1991",
                                          "Force,\nfrom 1991",
                                          "Overall,\nfrom 1991"),
                                 ordered=TRUE)

sub.bigfig$Variables.factor <- factor(sub.bigfig$Variables,
                                      levels=c("Common law",
                                               "Population (log)",
                                               "GDP cap (log)",
                                               "Polyarchy",
                                               "MID side B",
                                               "MID side A"
                                               ),
                                      ordered=TRUE)


sub.bigfig$Contrast1.factor <- factor(sub.bigfig$Contrast1,
                                     levels=c("Common law: Yes vs. No",
                                              "50M vs. 4M",
                                              "20M vs. 4M",
                                              "Population: 10M vs. 4M",
                                              "12K vs. 2K",
                                              "8K vs. 2K",
                                              "GDP cap (log): 4K vs. 2K",
                                              "8 vs. 2",
                                              "6 vs. 2",
                                              "Polyarchy: 4 vs. 2",
                                              "300M vs. 20M",
                                              "200M vs. 20M",
                                              "Milex cap: 100M vs. 20M",
                                              "3  vs. 0",
                                              "2  vs. 0",
                                              "MID side B: 1 vs. 0", 
                                              "3 vs. 0",
                                              "2 vs. 0",
                                              "MID side A: 1 vs. 0"),
                                     ordered=TRUE)


sub.bigfig$type1.factor <- factor(sub.bigfig$type1,
                                 levels=c("Overall",
                                          "Force",
                                          "No force"),
                                 ordered=TRUE)

sub.bigfig$type2.factor <- factor(sub.bigfig$type2,
                                 levels=c("from 1991",
                                          "from 1945"),
                                 ordered=TRUE)

cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

## knitr SUBSTANTIVE_WITH_POLYARCHY
figure5 <- ggplot(sub.bigfig, aes(color=Variables.factor)) + #
    geom_hline(yintercept=50, color=gray(1/2), lty=2) +
#    geom_linerange(aes(x=Variables, ymin=X2.5., ymax=X97.5.),
#                       lwd=1, position=position_dodge(width = 1/2)) +
    geom_pointrange(aes(x=Contrast1.factor, y=PI, ymin=X97.5, ymax=X2.5),
                    lwd=2, position=position_dodge(width = 3/5),
                    shape=21, fill="WHITE") +
#    scale_color_brewer(palette="Set1") +
    scale_color_manual(values=cbbPalette) +
    scale_y_continuous(limits=c(0,100),breaks=seq(0,100,10)) + #name="SATT",
    facet_grid(type2.factor~type1.factor) +
    coord_flip() +
    theme_bw() +
    #    ggtitle("") +
    theme(axis.text.x = element_text(size = 16,angle=45),#,
          axis.text.y = element_text(size = 16),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          legend.position="none",
          strip.text.x = element_text(size=20, face="bold"),
          strip.text.y = element_text(size=20, face="bold"),
          plot.title = element_text(lineheight=.8, face="bold",size=40)) +
    guides(color=guide_legend(nrow=2,byrow=TRUE,reverse=TRUE))

ggsave("gc-jg-project/kampala/paper/figures/figure5.png", plot = figure5,
    width = 14, height = 8, device = "png")


prob.mida.10.overall <- round(sub.bigfig[1,1])
prob.mida.20.overall <- round(sub.bigfig[2,1])
prob.mida.30.overall <- round(sub.bigfig[3,1])

prob.milex.10.overall <-  round(sub.bigfig[7,1])
prob.milex.20.overall <-  round(sub.bigfig[8,1])
prob.milex.30.overall <-  round(sub.bigfig[9,1])

prob.poly.24.overall <- round(sub.bigfig[10,1])
prob.poly.26.overall <- round(sub.bigfig[11,1])
prob.poly.28.overall <- round(sub.bigfig[12,1])

prob.gdp.24.overall <- round(sub.bigfig[13,1])
prob.gdp.26.overall <- round(sub.bigfig[14,1])
prob.gdp.212.overall <- round(sub.bigfig[15,1])

prob.claw.10.overall <- round(sub.bigfig[19,1])
